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SUMURY 




Research work on NASA Grant NAG 3-9 concentrated on development of 
numerical methods for the Euler equations and on development of analysis 
methods for these equations. Results of this work have been published as 
journal articles, AIAA papers, student theses, and MIT internal reports. 
The most important results were a streamtube Euler solver which combines 
high accuracy and good convergence rates with capabilities for inverse or 
direct mode solution modes and an analysis technique for finite difference 
models of hyperbolic partial difference equations. 

Two graduate students, Robert Bush and Michael Giles, were partially 
supported by this grant. Robert Bush received his S.M. degree in February 
1981 and his Ph.D. degree in September 1983. Michael Giles received his 
S.M. degree in September 1982. Michael was also supported for a summer 


visit to ICASE in 1983. 



TECHNICAL SUMMARY 


The earliest work on this grant was computational work on approximate 
factorization (AF) methods for Euler and Navier-Stokes equations and was 
performed by Robert Bush. The major new result from his work was recogni- 
tion that the optimum time step size for AF methods was not necessarily the 
largest stable time step. This result has now been substantiated by other 
analytical and computational results. A Jour, of Comp. Physics paper [1] 
and an S.M. thesis [2] were published during this phase. A copy of 
Reference 1 is included in the Appendix. 

The next phase of effort on the grant concentrated on a new analysis 
method applied to finite difference methods for hyperbolic equations, and 
was performed by Michael Giles. This work developed group velocity concepts 
for finite difference equations to explain spurious traveling wave 
solutions, dissipation and stability of inflow/outflow bound'ry conditions, 
and convergence rates. During this phase, a Jour, of Comp. Physics paper 
[3], an IMACS symposium paper [4], an ICASE contractor report [5] and a 
S.M. thesis [6] were published. References 3, 4, and 5 are included in the 
Appendix. 

The last phase of the grant work examined new solution methods for the 
Euler equations and was also conducted by Michael Giles. Major results 
from this work were a box-type method for the quasi-one-dimensional Euler 
equations, and a streamtube method for the two-dimensional Euler equations. 
Both methods are substantially faster than other comparable solution 
schemes. An AIAA paper [7] was published on the streamtube method, and an 
internal MIT report [8] was published on the box method. References 7 and 
8 are included in the Appendix. 
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Note 

Boundary Treatments for Implicit Solutions 
to Euler and Navier-Stokes Equations 


Introduction 

Implicit time marching schemes like those of Beam and Warning 1 1 ], Briley and 
MacDonald |2|, and MacCormack (1980) |3| generally have not been as robust as 
would be expected from a stability analysis for the pure initial value problem. 
Recently, Yee et al. [4 ) illustrated that a more general stability analysis, which 
includes the effect of boundary conditions, may explain some of the seemingly 
anomalous behavior of these schemes. T.ie major theoretical basis for this type of 
modal stability was established in a senes of papers by Kreiss 1 5. 6 1. Osher [7, 8 1, 
and Gustafsson et al. 1 9 J. 

Yee as well as Gustafsson and Oliger 1 1 0 1 considered the effect of inflow-outflow 
boundary condition formulations on the stability of a class of numerical schemes to 
solve the Euler equations in one space dimension. The characteristic feature of a 
subsonic inflow-outflow boundary is that a priori boundary values may be specified 
for only some problem variables, while remaining boundary values must be deter- 
mined as part of the solution process. Yee demonstrated a rather large disparity in 
stability bounds between the use of explicit or implicit extrapolation procedures and 
in general demonstrated that implicit extrapolation procedures had the least 
restrictive stability bounds. The intent here is to explore computationally the 
implication of this work for several two-dimensional Euler and Navier-Stokes 
simulations. 


Numerical Procedures 

The two-dimensional Navier-Stokes equations may be written in vector form as 


cU 8E cF cR cS 

+ — h-^ = -r— + ~T~- (I) 

cl cx cy cx cy 

The strong conservation law form may be retained under a general coordinate 
mapping as illustrated in Viviand (11]. All computations to be described were 
conducted in a mapped computational domain but for simplicity numerical and 
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Fig. I. Grid numbering schene for boundary condition formulat on. 

boundary condition procedures will be described in the simple two-dimensional 
geometry shown in Fig. I. 

A 1979 paper by Beam and Warning 1 1 2 1 outlined a solution scheme for systems 
of equations of form ( 1 ) which included most numerical schemes for which the modal 

boundary condition analysis has been conducted. This scheme uses the well- 

developed methods for ordinary differential equations as a guide to developing 
numerical methods for partial differential equations. The scheme presented combines 
linear multistep methods, local linearization, approximate factorization, and one leg 
methods. The shceme. a generalization of the scheme presented in 1 1 1. solves for a 
variable p(£> u which is equivalent to J u n in the class of schemes represented by the 
earlier paper. The earlier scheme is somewhat easier to understand as J u n is just the 
change in the solution from time level n to level n + 1. while p(E) u is a more general 
time differencing formula. 

The solution schemes chosen are implemented as 

(I + L m ,)Au* = RHS\ ( 2 ) 

(i + L' x ) AU n = JU*. (3) 

{/"** U n + AU\ (4) 

where RHS " is very nearly the finite difference approximation to the steady state 
equations, and L x and L y . are linearized difference operators representing a particular 
time and spatial differencing scheme. 

Full details of these operators are contained in |1|. If the spatial differencing is 
taken to be centered, the computational form of either Eq. (2) or (3) appears at each 
interior point as 

A 1AU1_ , + BlAL'l + C”AU% , = D”, ( 5 ) 

where A,, fl,, and C, are 4 x 4 matrices known at time level n, D • is the right-hand 
side vector at node point / known at time level n, and AU1 is the unknown vector at 
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node point /. The boundary points will be assumed to involve only the nearest two 
points in the x direction. 


A n 0 dU 0 + B n 0 AU t + C n 0 dU z = D n 0 . (6) 

The restriction to extrapolation along grid lines (actually transformed grid lines), is 
necessary to maintain 'he block tridiagonal form and avoids possible instabilities due 
to skewed extrapolation, see [13]. 

The full matrix equation will reduce to tridiagonal form if the first and nth 
equations are substituted into the second and the (n - i)th equations, for example. 

B\ = B X -A X A;'B 0 . (7) 


Boundary Treatments 


Inflow-outflow Boundary 

The finite difference algorithms studied usually require more boundary values than 
are required for the partial differential equations which they simulate. These extra 
numerical boundary conditions cannot be set arbitrarily and are usually determined 
through an extrapolation procedure. These extrapolation procedures may either be 
explicit, that is boundary values needed at a new time are determined uniquely from 
the old time level solution, or implicit, that is. the boundary values' are determined as 
part of the new time level solution. The analytical boundary conditions or the 
extrapolation quantities are usually not conservation variables but primitive variables, 
and a local linearization is usually required as part of defining the extrapolation 
procedure. 

Consider, for example, an implicit subsonic outflow boundary at which the local 
static pressure is specified as a boundary condition and all other variables are to be 
determined by extrapolation. Figure 1 shows a typical computational grid and defines 
the subscripts used. 


* 1 nn 

i-j r iJ 

IP) 

j p y* 1 

given. 

(8) 

= 2 l pu 

- 1 pu\ 

implicit space extrapolation. 

(9) 

i.i \ P v ) 

l.j-l \P v fiJ-2 




In order to complete the boundary formulation, all equations must be expressed in 
delta form and in terms of conservation variables. For the total internal energy this 
may be done through its definition 

£, = /»/(•/- \) + \(pu) 1 /p + (v) i /p. (10) 

Since the relations between conservation variables are nonlinear, some linearization 
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step will be necessary before the boundary condition formulation may be used. We 
choose to introduce our linearization step here as 

J£ t = (£?*'-£?) = (l/(y- + + 

+ (JuJf.Ju 2 . Jv 2 , ApAu, ApAv). (11) 


If terms of order Judi- are neglected, the error is equivalent to the linearization error 
of the interior point scheme. We may express the transformation from boundary 
variables to conservation variables as 


U'.j 




I 0 0 

0 1 0 

0 0 1 

*(u 2 + r 2 )" u r v" 


0 

0 

0 


1/(7- 1) 





( 12 ) 


we shall in general denote transformation from conservative to primitive variables as 


= T'jdUjj. 

The extrapolation conditions for W Lj are 



(13) 


(14) 


or 




(15) 


The final equations relating the boundary conservation variables and the interior 
conservation variables are 




or 


(15) 


i.fii j_ i^Wu - 1 + ( 1^1 

With the definition of and given in Eq. (15). T, J _ i and T u _ 2 are identity 
matrices. 
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An explicit outflow boundary treatment was constructed using 

pn *l_pn given, 

Ipu ) Mu) . (18) 

\PvJu \pvlu-i 


and setting G u _ , = H u _ 2 = 0. 

In forming Eq. (9), we choose to extrapolate the local momentum flux rather than 
a specific primitive or characteristic variable; choice of other extrapolation variables 
would alter only the transformation matrix T. y . Extrapolation of the momentum flux 
is somewhat arbitrary, but its choice did not affect the accuracy of the computational 
results to be presented. 

Solid Wall Boundary Procedures 

The boundary treatment procedure illustrated for inflow-outflow boundary are 
easily extended to cover solid walls in either inviscid or viscous flow situations. Here 


or 




Ap \ 

/ y/T 


P/T 

Apu \ i 

1 yu/T 


Pu/T 

Ape U 

yv/T 


pi IT 

AEy ) ’ 

1 1 + 1 

yr 

1 pq 

\ /’ — 1 2 

T 

2 T 


AU 0J = 

V" 

‘V/ 

AW 0J . 


Tlr P <? 



( 10 ) 


( 20 ) 


where q is the velocity parallel to the wall and S is the wall slope. For the inviscid 
flow examples cP/c\\ cT/cy. and cq/cy are set equal to zero, while, for the viscous 
flow examples v. u, and cT/cY are set iqual to zero and cP/cy is equal to 
4/3yu(r J /0' : )(i')- All derivatives are evaluated by one-sided finite difference formulas. 

As indicated by Buggeln et al. [14), an ADI type procedure requires boundary 
conditions for the intermediate step. Usually the intermediate step was in the y 
direction and the boundary conditions were applied as if the intermediate results were 
physical quantities, that is. the boundary conditions of Eq. (19) were applied to the 
quantities AU * of Eq. (2). 

Explicit wall boundary treatments are generated by applying the primitive variable 
form of Eq. (19) and forcing the correction matrices to be zero. 


Numerical Results 

Three geometries were selected for detailed study: an inviscid supersonic diffuser 
with w'eak oblique shock, supersonic in/supersonic out; an inviscid supersonic 


BOUNDARY TREATMENTS 


307 



Fig. 2. Computational grid for weak shock difTuser calculations. 

diffuser with a strong normal shock, supersonic in/subsonic out. and a viscous super- 
sonic diffuser with weak oblique shock illustrating a shock-boundary layer 
interaction. Sketches of the geometries are shown in Figs. 2-4. Solutions Tor each 
geometry were run to steady state fcr a range of time step sizes. For convenience, 
time step sizes are reported in terms of x and y CFL numbers 

(CFL,) = max(J/(u + c),JA. r fJ ). (21) 

(CFL), = max(Jf(r + c\JAy,J. (22) 

The time step size was uniform over each calculation which results in nonuniform 
CFL, and CFL, numbers. The maximum value of each is reported. Sample 
convergence history plots are shown in Fig. 5 which shows the log of the value of the 
point maximum steady state residual 


. SSR = cEjcx - cR/cx + cFjcy - cS/cy (23 ) 

plotted against the iteration number. A solution was not termed stable unless the 
residual converged to the machine accuracy, about 1 x 10'*. All calculations used a 
32 bit floating point word size. 

Each geometry calculation was run with fully explicit extrapolations. Au = 0. and 
with fully implicit extrapolations: the results are summarized in Table I. The most 
interesting of these results are shown in Fig. 5. At a tipe step size corresponding to a 
CFL, number of 15. convergence was rapid and very nearly monotonic in time. At 
smaller time step sizes, the convergence was slower but nearly monotonic. At a CFL, 
of 45. convergence rates initially appeared to be faster than for a CFL, of 15. but the 
final residual values oscillated sugnificantly about its minimum value. At a CFL, of 
90. the convergence rate was substantially slower than at a CFL, of 15. and at larger 
CFL, values the solution diverged. 

The results for 'he strong shock diffuser can reasonably be compared to those of 
Yee et al. [4|. They reported a CFL number stability limit between 10 and 20. while 



Fig. 3. Computational grid for shock-boundary layer calculations. 
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ITERATION l 10 


Fig. 5. Convergence history for strong shock diffuser calcuimon. 
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we found stability limits between 90 and 130. Thus the analysis in one space 
dimension does appear to provide a sufficient condition for stability, but it may not 
provide a close approximation to the stability limit. It is essential, however, to 
emphasize that the largest convergence rates were observed at time steps 
corresponding to CFL numbers of order 10 and that only a marginal computational 
time advantage for the implicit boundary formulations was observed. 

The results for the shock-boundary layer calculation are very interesting but they 
demonstrate a substantial computational advantage for the implicit solid wall 
conditions, not for the inflow-outflow extrapolation. Here the stability boundary and 
the best convergence rates were observed at time step sizes corresponding to CFL, 
numbers of 5 to lO. When using the ir solicit wall conditions, the algorithm stability 
appeared to be independent of grid spacing in the normal direction as might be 
hoped. When using the explicit wall condition, the algorithm stability was limited to a 
CFL,. number of about 500. 


Conclusions and Discussion 

While it is difficult to generalize from only a few test camples, it .s apparent that 
a better appreciation of the role that boundary treatments play in implicit algorithms 
has a'lowed the development of far more robust Beam and Warming type solvers. For 
both explicit and implicit boundary treatments, we were able to compute solutions 
accurately with time steps 50 to 100 times ger than explicit time limits while 
retaining the ability to choose rather arbitrary initial conditions. In many cases, our 
limiting time steps for the two-dimensional test problems were in fact larger than the 
limit which a one-dimensional analysis would suggest. 

The most important computational result we observed was that while an improved 
appreciation of boundary treatments did allow very large time step sizes to be used, 
the largest convergence rates to steady state were observed at relatively small time 
step sizes. For the two-dimensional test problems, the best CFL, numbers were of 
order 10. not of order 100. One-dimensional test examples showed no such 
converg nce rate behavior. Present!, unpublished analysis by Abarbanel el ul. |I5| 
has linked this behavior to the approximate factorization form of Eqs. (2) and (3). 
This efTect now seems to be setting the time step sizes for our viscous How 
computations and new work should focus on methods 'or overcoming this limitation. 
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Abstract 




An asymptotic approach is us«d to analyze the propagation and dissi- 
pation of wavelike solutions to finite difference equations. It is 
shown that to first order the amplitude of a wave is convected at the 
local group velocity and varies in magnitude if the coefficients of the 
finite difference equation vary. Asymptotic boundary conditions coup- 
ling the amplitudes of different wave solutions are also derived. 
Equations are derived for the motion of wavepackets and their interac- 
tion at boundaries. Comparison with numerical experiments demonstrates 
the success and limitations of the asymptotic approach. Finally an 
asymptotic global stability analysis is developed. 


Notation 


5 a . 

X 


a. - u. 

1 +' 3 


u 0 . 

X 3 + ! 


i ‘ V 


V 


A U - U 4 - IT 
x J 1+1 j 


7 

x 


u. - u. - u 
J J 


j-1 


E U. - U. 
mx j j+a 


When there are several independent variables the subscript on the 
finite operator denotes the direction of the shift, differencing or 
averaging. For example, 


if 


U - u ( x , t ) 
J j n 


then 


„n 

6 x V* 


u n - u n 
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A » complex conjugate of A 

Re (A) - Real component of A 
Im(A) • Imaginary component of A 


I. Introduction 

Methods for analysing dispersive partial differential equations are 
veil established. Using Fourier decomposition and asymptotic evaluation 
of integrals, or by direct asymptotic expansion, [5,9] it can be shown 
that the energy propagates at the local group velocity. Ray theory 
[5,9] then treats wavepackets, localized wavelike disturbances, as 
particles and derives 3imple o.d.e.'s for their motion. This paper 
applies the techniques to the analysis of numerical wave propagation in 
finite difference equations. Due to the discretization the numerical 
waves are always dispersive even if the analytic system being modeled is 
nondispersive. Until recently the importance of the group velocity in 
analyzing finite difference solutions does not seem to have been 
recognized. Kentzer [4] discusses the role of group velocity and shows 
'that in many common schemes the numerical group velocity at high 
wavenumbers 1s in the opposite direction to the analytic group velocity. 
Vlchnevetsky and Bowles [8] derive reflection coefficients for the 
Interaction of waves at boundaries, and present several illustrative 
numerical examples. Trefethen [6] provides a group velocity interpreta- 
tion of the stability theory of Gu3taffson, Kreiss and Sundstrom [3], 
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and in a forthcoming paper [7] will derive rigorous conditions for the 
P-stability [10] of two-boundary problems. In the stability analysis in 
this paper we use P-stability, which is concerned with stability in the 
limit t«->, rather than GKS-stabilty which is concerned with stability in 
the limit At-0. Reference [2] contains further details and numerical 
examples of the work in this paper. It also includes' a more general 
global stability analysis which allows for variable coefficients in the 
finite difference equations, and in the case of constant coefficients 
reduces to the exact stability analysis of Beam, Warming and Yee [1]. 

The approach we use is an asymptotic one in which a wave solution is 
expressed as a product of a complex amplitude and an oscillatory phase 
function whose frequency and wavenumber may also be complex. The asymp- 
totic assumption, or approximation, is that the length scale for varia- 
tions in the amplitude and wavenumber is large relative to a mesh cell 
length. An asymptotic expansion leads to a local dispersion relation 
relating the wavenumber to the frequency. The first order ten.s produce 
an equation for the amplitude in which the local group velocity appears 
as the velocity of convection of the amplitude. AI 30 there is a 
t variation in the magnitude of the amplitude if the coefficients of the 
finite difference equation vary. All of the wave solutions with a given 
frequency and different wavenumbers are coupled at the boundaries by 
asymptotic boundary conditions. If there are only two waves per 
frequency then this reduces to the amplitude reflection coefficients 
computed by Vichnevetsky and Bowles [8]. 
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The following section develops a theory for the notion of wave- 
packets which are wavelike disturbances of finite length and constant 
frequency. Using the techniques of classical ray theory [5,9], these 
can be treated as particles and sinple o.d.e.'s can be derived to 
describe their notion and the change in their energy. When they reach 
the boundary they are reflected into wavepackets of a different wave- 
number but the same frequency and the energy of the reflected wavepacket 
can be calculated from boundary reflection coefficients. The last 
section derives a global stability analysis in which the usual Fourier 
stability analysis is modified to calculate the effects of non-periodic 
boundary conditions and slowly varying coefficients. This analysis is 
then used to calculate the spectral radius of the backward Euler method. 


II. Asymptotic Amplitude Analysis 
Asymptotic Amplitude Equation 

Consider a general linear homogeneous finite difference equation 
with variable coefficients, 

u" - 0 (1) 

where 

L = 7" C (x) E E (2) 

J mp mx pt 

m,p 
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If the coefficients C ere constants then 

Dp 


• exp[i( j>-n3) ] 


( 3 ) 


is an exact solution of ( 1 ) provided 




n,p 


C exp(i (m«-p^ ) ] ■ 0 
op 


(4) 


This equation is called the dispersion relation between wavenumber $ 
and frequency 1 . If the coefficients are not constant then U can be 
expressed as, 


• A(j,n) expliM j,n) ] 


(5) 


where A(j,n) is a slowly varying amplitude and y(j,n) is the phase of 
the wave and is related to the frequency 1 and wavenumber ? by 


3n 




3 ? 

— * t 

33 


( 6 a, b) 


The asymptotic approximation which is made is that the length scale 
I>A and time scale Ta for variations in A and the length scale for 
variations in 9 are ouch greater than 1. Substituting (5) into (2) and 
expanding A and y in Taylor seies about a point (j,n) yields, 


r— » 

L. u" - exptly] / C (j) exp(i(D9-pi7 ) ] 
j J < 


m,p 


, 3 A 3 A im 2 34 

33 P 3n 2 3 j 


♦ 0< AL “2, AT " 2 , AL " 2 ) 
A A ^ 


(7) 


To satisfy equation (1) the amplitude A(j,n) must satisfy , 
a 0 U ,fl, j) A ♦ a ; U,2i j) U.3, j) ~ ♦ a, (i ,.3. j) A - 0 (8) 
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where , 


*, (* .3,3) 

*1 (<*> .3,3) 

«i (« i 3 . J ) 


V C (3) exp[l(m®-p:,3) ] 

L , Dp 

n,p 


iii’l . 

j 1 i 3 


const 


f4i.' 

[ot n,j const 


»,(«,3,3) - - - 7:777 


i f b 2 a % 

2 U® 2 ; 3,3 const 


(9a) 


(9b) 


(9c) 


(9d) 


Because of the asymptotic assumptions (8) can only be satisfied if 


*,(5.3.3) - 0 ♦ 0( lj- , .l A - , i T A - 1 ) 


( 10 ) 


This is the asymptotic form of the dispersion relation between $ and 
Q and will usually be satisfied by setting ao identically equal to zero. 
$ is now a slowly varying function of 3 due to the slow variation in the 
coefficients. Neglecting the second order terms and dividing by ai gives 
the asymptotic amplitude equation. 


3 A 3 A 

— ♦ r a — - ■ c A 

3n 9 3] 


, where r g « a, / a t 


and 


- (a 


35 


~ ♦ a ) / a 

5 33 J ' 1 


( 11 ) 


(12,13) 


Differentiating Eq.(33) with 3 held constant gives, 


da, 


(3a.' 

■ 


'3*,' 


*a. 1 ia i 

[33 <5,3 const 1 3$ 3,3 const 


d® - 0 


(U) 


Hence , 
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Bt, j const • *i ' *> ’ r , <,5> 

Thus ths amplitude A of the wav* is convectad at tha local group 
velocity. 

Asymptotic Boundary Conditions 

Tha general solution of Eq.(1) is a sum of waves with different 
constant frequencies 3 and slowly varying wavenumber $ and amplitude A. 

— J* j 

u" • ■ ) A (j.n) exp[i( f * { r ) dr, -n3 ) ] (16) 

j / , ° J o 

3 m-1 0 

The outer su mm ation is over different values of 3, and the inner 
summation is over the M different values of $ which satisfy the 
dispersion relation for each 3. For each 3, a the amplitude A satisfies 
the asymptotic amplitude equation on the interior of the computational 
domain Independent of all the other waves. All the waves of each 
frequency are however coupled by boundary conditions. 

Suppose one of the finite difference boundary conditions at j-J is 

B u" s ) D. E, E u" • 0 (17) 

J L-j lp lx pt J 
l.P 

Performing the same expansion as in the derivation of the asymptotic 
amplitude equation, retaining only the leading terms, and equating the 
coefficients of exp(-ift) for each 3, the boundary condition becomes, 
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( 18 ) 


M J 

r b(fi,$ ) A (J) axp[i ^ * dr ] ■ 0 

Z | OB j O 

n- 1 0 


where, 


b(n,$ a ) - ^ d •xpddi-pn)) 

i.p 


( 19 ) 


There art similar asymptotic boundary conditions at j-0. 


III. Rav Theorv And WaveDaeket-Particles 


In addition to tha asymptotic approximations mad* aarliar this 
saction assumes that for all real wavenumbers % , the frequency 3 is real 
for all j and so the group velocity r g i 3 real. 

A Lagrangian-type total time derivative is defined by, 


d_ 

dn 




( 20 ) 


so 


- r 

dn 9 


( 21 ) 


dA 

dn 


c A , 


( 22 ) 


• nd 

A general initial value problem for a wave of frequency 3 and wave- 
number $(3>J) can be solved by integrating these equations from given 
initial conditions. 

A wave packet is a wave for which tha amplitude A is non-zero on only 


Page 10 



a small part of tha domain. Tha energy la defined to be, 


E(n) 


•/ 


A(x,t )| dx 
n 


J 

/ 


| A( j , n) | •gy d j 


( 24 ) 


Differentiating this definition, and using (22), yields, 
dE . - fdx'-’a , dx. . _ 

3T ■ 1 c * c tj (r 9 dj> 1 E 


(25) 


Thus equations (21) and (25) describe the motion of a wavepacket 
particle in the interior of the computational domain. When the wave- 
packet reaches a boundary it is reflected as one or more wavepackats 
with the same frequency but different wavenumber. For the case in which 
there are just two wavenumbers corresponding to tha same frequency the 
ratio of the reflected energy E, to the Incident energy is given by, 


h 


fq ( * . , ) 

*g(s i J 


where the amplitude reflection coefficient R is defined by, 


R. ■ 


A j ( J , n * 


J A^J.n) 

and is determined from the asymptotic boundary condition. 


(26) 


(27) 


Example 


The example is the solution of the modal convective equation, 


at ax 


( 28 ) 
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using a trapezoidal method, 


* T V» 


“t 3 x 


1 

2 r H 


u 7 

t X 


)■? 


n *i - 


(29) 


with r . - — - C — — (30) 

Vi’*! 

Figure 1 ahowa the solution U(x,t) corresponding to a uniform grid 
0<j<2C0 with r*1, and initial conditions corresponding to a wavepadeet 
approximately 20 mesh units wide. Comparison of the heights of the 
wavecrests a-e at time levels 60 and 120 shows that the phase velocity, 
the velocity of the wavecrests, is greater than the group velocity, the 
velocity of the wavepadeet. 

Figures 2 and 3 show comparisons of the wavepacXet theory w.th 
numerical experiments. Zn each case "experimental" values fer X(n), the 
position of the wavepadeet, and E(n), its energy, are obtained by 
solving the finite difference equations and "predicted" values are 
calculated by solving the wavepadeet equations. The Initial wavepadeet 
in each case is similar to that in the previous example. 


Zn the first case r varies exponentially from 0.0S at j-0 to 0.2 at 
f j»200 and 3 * 0.04 , Figure 2 shows X(n) and ln[E(n)] both predicted and 
experimental. This example shows the movement of a wavepadeet and the 
change in its energy due to the variation in r. The agreement between 
the predicted and experimental values is excellent. The energy of the 
analytic solution is constant so the wavepadeet theory has successfully 
predicted almost all of the change in the numerical energy due to the 
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nonunifora grid 


Zn the second csss r is constant and eq. v to 1.0 and 3-0.3 . 

Figure 3 shows X(n) and ln(E(n)]. This example illustrates the effect 
of the downstreaa boundary reflecting e wevepacket with reduced energy. 
Because of the finite length of the wevepacket the drop in energy la 
smeared and X(n) does not quite reach 1.0 . Again the agreement is 
excellent with the wevepacket theory accurately predicting the energy of 
the reflected wavepacket. 


IV. Asymptotic Stability end Convergence Analysis 

Zn this section it is assumed that there are two wavenumbers 
corresponding to each frequency, and that if one is real then so too 
is the other. Examples of method* satisfying these conditions are the 
trapezoidal method applied to the a>xiel convectlvn probl# i with variable 
CFL number r, and the whole class of Beam-Warming schemes applied to the 
model convective problsm with constant CFL number. 

The normal Fourier analysis assumes constant coefficients and 

r 

periodic boundary conditions and derives eigenf requencias 3 ( 3 ) where 3 
is a real wavenumber satisfying the periodic boundary conditions and 
Q ( a ) la the corresponding frequency given by the dispersion relation. 

The common use of Fourlnr analysis to predict the stability of problems 
with nonperiodic boundary conditions implicitly assumes that 3 ( 3 ) is a 
close approximation to the true eigenf requency. This section follows 
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that assump cion , calculates a correction 3* to the Fourier frequency 
aU) due to the boundary conditions, and then determines the validity of 
the assumption based on the asymptotic errors. 


r i 

u” - A 2 ( j ,n) exp i j d£ - 

J l 0 


( 1 


* A 2 ( j , n) exp 


j ♦, dr -inn (3 


1) 


If the true eigenf requency is n+fl* then A«exp(-inn') and so, 

A (j,n) » exp(-ina') A ( j , 0 ) , m-1,2 (32) 

n m 

Substituting these expressions into the asymptotic amplitude 
equations to evaluate the time derivative and then integrating the 
resultant o.d.e. gives, 


A (J,0) 
m 


A (0,0) exp 


U r g - 


dj , m-1 ,2 


(33) 


The two asymptotic boundary conditions then become in matrix form 

I A, (0,0) I 


(34) 


A 2 (0,0) i 


A non-trivial solution exists only if det^B) ■ 0 and this leads to 
'the following equation for 3*. 


exp 


' ? 


i 


in* j [r g (# l , j) ] -1 -[r g U 2 , j )]” 1 dj j - 
0 1 


b 2 ( ft i £ ^ iii j ) 


w 

1 / *a’’i * xp 


j ■- - * dj 

0 ‘ J * 7 ‘ j 


(35) 
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If ■>, are chosen so that the r.h.3. is real and positive then 


n ' 



b 2 ( a,» : ) b ; (a ,e l ) 

bj ( .1 , $ ^ ) b^ (.1 , i j ) 




where, 

J 

N - f C r g ( * t , j ) ] “ 1 - [rgUj.j )]" 1 dj 
0 


(36) 


(37) 


Thus the frequency a resulting from a normal Fourier analysis is 
corrected by an amount a ' due to boundary conditions and variable 
coefficients. This approach, using a as an initial approximation to the 
actual eigenfrequency , is valid provided the asymptotic error is small. 

The asymptotic error is 0 (L a “ 2 ,T a “ 2 ) ■ 0(J" 2 ,a'^) so provided r^<<J 
N >> I and hence a' << 1 except near frequencies for which 

bj (a ,s 7 ) b. (a ) 

bj (a.ij ) b t (a.Oj ) 


is zero, or infinite, which usually occurs at n»0. However these 
frequencies are heavily damped by the boundary conditions and so an 
accurate estimate of their eigenf requencies is not essential. This 

P 

method gives accurate asymptotic values near the critical frequencies 
which are least damped and which therefore determine the overall 
spectral radius of the scheme. 
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Examole 


This example is the backward Euler method applied to the model 
convective problem with constant CEL number r and space extrapolation at 
the downstream boundary. The finite difference equation is, 

( 7 t ♦ r “x , x) U "*' 1 ' 0 ' 1381 

and the dispersion relation is, 

a 9 - 1 - exp(i.l) * ir sin(o) - 0 (39) 


After carrying out the calculations the frequency correction 3* is 
found to be [2], 


Q* - - 


ir cos ( ? ) (1 - ir sin( : ) ) 
2J ( 1 ♦ r i sin i (« ) ) 


log[cot( ;/2 ) ] 


(40) 


Thus the effect of the boundary conditions i3 to greatly accelerate 
convergence at low wavenumbers while having little effect on the higher 
wavenumbers. The spectral radius X is 


\ - max \ exp( -i(fl*3 * ) ) 
ft 

r log(J) 


1 - 


4 J 


for J>>r 


(41) 


r 


V. Conclusions 

The validity of the asymptotic approach developed in this paper is 
demonstrated by the numerical results in sectic III. The limitations 
of the wavepacket theory are due to the asymptotic approximations 
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involved in treating the wavepacket as a particle. The stability 
analysis in section IV uses fever approximations and so the asymptotic 
errors will be substantially smaller. 

The calculation of the asymptotic amplitude equation and asymptotic 
boundary conditions for a particular case is no more difficult than a 
normal Fourier analysis. For applicable cases the vavepacket theory and 
the stability analysis are straightforward. In more complex cases the 
main benefit from the theory is the Insight given by the asymptotic 
amplitude equation and boundary conditions. The amplitude equation 
gives the group velocity and the effect of varying coefficients which is 
of great interest since in 2-D cascade geometries cell lengths can vary 
by factors of up to 100 in inviscid calculations and 1000 in viscous 
calculations. The asymptotic boundary conditions give the amplitude 
reflection coefficients which provide a practical criterion for choosing 
the best numerical boundary conditions. 


Acknowledanent3 


This research was supported by NASA Lewis Research Center under 
p grant No. NAG3-9 with Technical Monitor Or. R. V. China, and by a 
Scholarship for M. Giles from the Kennedy Memorial Trust. 


References 

1 R. Beam, R. Warming, H. Yee, "Stability Analysis for Numerical 
Boundary Conditions and Implicit Difference Approximations of 
Hyperbolic Equations", Proc. NASA Symp. on Numerical Boundary 
Procedures, 1981, pp. 99-207. 


Page 17 



2 M. Giles and W.T. Thompkins Jr. , "Asymptotic Analysis of Numerical 
Wave Propagation In Finite Difference Equations" , MIT Gas Turbine 
and Plasms Dynamics Laboratory Report No. 171 , February 1983 

3 B. Gustafsson, H-O. Kriess, A. Sundstrom, "Stability Theory of 
Difference Approximations for Initial Boundary Value Problems II", 
Math. Como. 26 (1972), pp. 649-686. 

4 C.P . Kentzer , "Group Velocity And Propagation Of Numerical Errors" , 
AIAA Paper No. 72-153 

« 

5 J. Liqhthill , Waves In Fluids , pp 237-260 , Cambridge Univerity 
Press 1978. 

6 L.N. Trefethen, "Group Velocity Interpretation of the Stability 
Theory of Gustafsaon, Kreiss and Sundstrom", Jour. Comp. Phys. 49 
(1983), pp. 199-217 

7 L.N. Trefethen, "Stability of Finite Difference Models Containing Two 
or More Boundaries", to appear 

8 R. Vichnevetsky and J. Bowles , Fourier Analysis of Numerical 
Approximations of Hyperbolic Equations , SIAM Studies In Applied 
Mathematics 1982. 

9 G.B. Whitham , Linear And Nonlinear Waves , chapter 11 , John Wiley & 
Sons 1974. 

10 H. Yee, R. Beam, R. Warming, "Stable Boundary Approximations for a 
Class of Implicit Schemes for the One-dimensional Inviscid Equations 
of Gas Dynamics", AIAA-31-1009-CF , AIAA Computational Fluid Dynamics 
Conference, Palo Alto, June 22-23, 1981 





In (E) 


NUMERICAL 

EXPERIMENT 

WAVEPACKET 

THEORY 



oo 



Figure 2. Position and Energy of Vfavepacket: Effect of 

Non-Constant CTL Nunber 





In (E) 


NUMERICAL 

EXPERIMENT 


WAVEPACKET 

THEORY 




Figure J. Position and Energy of Wavepacket: Reflections 

at Boundaries 



INTERNAL REFLECTION DUE TO 
A NONUNIFORM GRID 



ORlGiMA- I '* - 

OF POOR QUW-llY 


Nicna«-i 3»l«a 
W. T. Theme* Jr. 

Massacnusatts Institute it Teehnoloqy 
Cambridge. IA 02139 
Tel: 617-253-2276 


Abstract 


•■•■pita of wave- trapplnq , tho Internal 
raflactlon of nuaarlcal wavaa duo aolaly to 
oar la tlona In qrld atrotehinq. In both 
aaanplaa tha analytic aquation la tha acalar 
convactlon aquation which la non-dl apo r a 1 va 
and for which aach rourlar conponont of a 
general dlaturbanca propaqataa at tho aano 
oaioclty c>0. Howovar, nuaarlcal approxina- 
tlona of thla aquation ualnq 3-polnt apatlal 
dlf f aranclnq and trapoaoldal (or Crank- 
Nlcholaon) tlna Intaqratlon hava tha proparty 
that for a qloan fraquancy thara la a ainuaun 
raqulrad apatlal raaolutlon ax for travallnq 
wava aolutlona. In both tha aaanplaa praaan- 
tad, an Initially wall raaolvad wavapacaot 
propaqataa towarda a raqlon of tha qrld In 
whicn tha apaclnq n la lncraaolnq, until It 
raachaa tha point at which tha apatlal raaol- 
utlon raachaa tha critical value and ita 
q roup velocity la taro. It la than raflactad 
and baconaa a wavapacket with wavalenqth laaa 
than 4 node polnta travallnq with naqatlvo 
qroup velocity. It travala throuqh tha well- 
raaolvad raqlon until once aqaln It raachaa a 
raqlon of Inadequate qrld raaolutlon, and la 
than raflactad back Into a wavopackat with 
wavalenqth qraatar than 4 node polnta and a 
poeltlve qroup velocity. Tha difference 
botvoen tho two aaanplaa llee In tha detalla 
of tha apatlal dlf f aranclnq, which cauaaa no 
qualitative chanqe but qraatly affaeta tha 
•anerqy* of tha wavopackat durlnq tho 
oaclllatlon. 

Thaaa two caaaa are analyied ualnq a 
pravloualy darlvad aaynptotlc analvata which 
calculatea alnpla o.d.a.'a for tha notion and 
tha anorqy ehanqa of wavopacketa travallnq 
throuqh nonuniforn donalna. Caaplta tho 
praaanca of tha turnlnq polnta at which tha 
alnpla aaynptotlc analyali la not etrlctly 
valid, tha aqraonont batwaan tha nuaarlcal 
experlejnta and tha thaorotlcal analyala la 
excellent. 
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I . Introduction 


aany nunarical ca leu la 1 1 on a . Por oiaaple, In 
calculatlnq two-d 1 nana 1 ona 1 tranaonlc flow 
ovar an airfoil, tha qrld apaclnq ax win 
uaually bo vary anall nacr tha laadinq adqo 
of tha airfoil to raaolvo larqa q r all a n t a , 
and ardund tha ahock to Halt tha arrora dua 
to nunarical dlaalpatlon, while In tha far 
flald ix oftan baconaa vary larqa. Thara la 
nuaarlcal evidence that thaaa nonuniforn 
qrlde can cauaa aoaa problana. Croach and 
Orataq (21 found apurloua non-phyaieal 
Intarnal raflactlona dua to qrld atratchlnq. 
Ranee thara la lntoraat In analvtlnq alnpla 
nodal problana to qaln Inalqht into tha 
dlfflcultlaa. 

Thara are two llaltlnq caaaa of non- 
unifora qrlda. In tho first there are two 
unlforn qrlda with different spatial raaolu- 
tlon si Joined by an interface. This ease 
has bean analyttd by Vlchnavotaky [6] and 
Trafathan [51. To aunnarlia thalr flndlnqs. 

In qenaral a wave Incident on tha Interface 
produces a raflactad wava In addition to a 
trananittad wave. If tha wavolanqtn of tha 
wave la auch larqar than n on both aides, 
than tha raflactad wava haa a vary snail 
anplltuda. If tho wavalenqth la of tha aano 
order aa n than up to loot of tha nunarical 
enarqy la raflactad fron tha lntarfaea. 

In tho second caaa ax varlaa slowly ovar 
a lanqth scale auch larqar than ax. It la 
this case which la conaldarad hors u.inq a 
pravloualy derived aaynptotlc analyala ( 1 ) . 
This analysis la also applicable to finite 
dlfferanca aquations In which tha non-unlfora 
coefficients are duo not to non-unlfora qrlda 
but to elowly varylnq analytic aquations. 
First ve review tha theory for the aaynptotlc 
analysis and than wa praaant two nunarical 
aolutlona of tha nodal convection aquation 
ualnq 3-polnt spatial dlffarenclnq and trape- 
zoidal (or Crank-Nicholaon ) tins Intaqratlon. 
■acauso of tha particular choice of variable 
qrld apaclnq ax and tha fraquancy of tha 
wavopackat ehoaon for Initial conditions, tha 
exaaples both axniblt wava-trapplnq in which 
tha wavapackat'a notion la confined to a 
central portion of tho donaln and Its posi- 
tion, wavanunoar and anorqy qo throuqh a 
periodic oaclllatlon. It la shown that tha 
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aayaptotlc theory iceurittiy predict! this 
behavior. 


I ! . k a v 1 a w of Aaviantotlc >nilv»li 


Conildar a qenaral Umar, hoiodimoui, 
(Inlta difference aquation with varlaola 
coa (flclanta, 


U* ■ 0 < l ) 

vha ra 

l i * Ts c -p ( j ’ m 

■.p 

Pollowlnq tlia aayaptotlc approach of 
claaaleal analyala of wava propaqatlon In 
nonunlfori dleparelve aadluaa, wa conildar 
tha trial aolutlon. 


tl" ■ A( J , n ) aap ( IT ( J , n I | 


( 3 ) 


whara A(),n) la a alowly varylnq aaplltuda 
and f<3,n) la tha phaaa of tha oawa and la 
ralatad to tha fraquancy a and wavanuaoar a 
by 


21 . 

in 


-a 



(4.3) 


Tha aayaptotlc approalaatlon which la 
aada la that tha lanqth acala L for varla- 
tlona In A and a la auch qraatar than 1. 
Subatltutlnq (3) Into (1), eapandlnq A and f 
In a Taylor aarlaa about a point l),nl, and 
naqlactlnq taraa of ordar l*' , producaa. 


whara , 


*» A * *'71 * * 


’•IT * *■ 


i : 
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L c .p ,3) 

■ . p 

aap ( 1 ( aa -pQ ) 1 

(? ) 


-1 —9 

1 1 ' » 3 

. a . • -- ■ . J . 

( 8- 10 ) 

10 * * 

la 

• 2 1 a 1 


All of tha taraa In (6) aicapt for tha 
flrat ara of ordar L"', and ao wa raqulra 
that tha O(t) tara In a 4 la Identically equal 
to zero. Since 1 la aaauaad to be conatant, 
bacauaa tha coafflelanta In (2) ara not 
tlae-varylnq , thla condition daflnee tha 


dlaparalon aquation a s a(0,)). 

equation 

( 6 ) 

can than ba raarranqad to fora 

tha aayaptotlc 

aaplltuda aquation. 
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r * " conat 


(12) 


-• /• 

o 1 

(13) 


Thua tha aaplltuda A of tha wave 1* 
convactad at tha local qroup velocity and 
qrowa or daeaya In nonunifora caaaa in which 

tha wawanuabar a la a function of ). 


Ha now raatrlct our attantlon to non* 
dlealpatlva nuaarlcal achaaaa for aodalinq 
hyperbolic lyataai. In thaaa caaaa tha jroup 
velocity for tha nuaarlcal achaaa la real, 
and It la convenient to define a Laqranqlan- 
typa total tlaa derivative, 


d 1 

dn 1 in r 9 } J 


( 14 ) 

ao that, aa In claaaleal ray theory 

(•) . 

, a 

aat of aquatlona aay ba written for 
propaqatlon alonq raya. 

t ho 

wava 

Si " r 9 • 


1 13 ) 

»-■* • 


(16) 

H - r„ 1± . 

dn 1 i J 


(17) 

A qanaral Initial valua problaa 

for 

a 


wava cf fraquancy Q and wavanuabar #(0.3) can 
ba aolvad by inteqratlnq thaaa equatlone troa 
divan Initial condltlona. In ( 4 1 Trafathan 
derived tha klneaatle ray aquatlona (HI. (IT) 
froa tha dlaparilon relation for an anlao* 
tropic 2*0 caaa In which tha qrld waa unlfora 
but tha analytic coafflelanta wera not. Coa- 
putatlonal aiparlaanta conflraad tha predlc- 
tlona of the ray theory. In thla papar wa 
are lntaraatad In tha notion of wavepackata, 
which ara wavaa for which A la taro on all 
but a aaall part of tha doaaln. Tha anerqy 
la daflned to be, 

*J 

(In) • j | A ( * , t ^ ) | 2 da 

“o 

J 

- | | A ( J , n ) | 2 jy dj Ml) 

0 

Olf ferenilatlnq thla definition, and 
ualnq (16), y 1 a Ida . 

Si • ( 2 *' ,t, *(ST)’'h ,r ^” E n ” 

Thua aquatlona (IS), (IT) and (19) 
daacrlba tha notion of a wavapackat In tha 
Interior of the eoapu ta t Iona 1 doaaln. 


III. Eaaaplaa and Analyala 


•oth of tha nuaarlcal eaaaplaa ara 
•olutlone to tha nodal convective aquation 
with conatant poaltlve velocity c. 

Ju 


e ■“ - 0 

It > a 


( 20 ) 


■oth achaaaa uaa trapezoidal tlaa lntaqr- 
at Ion and 3-point apatlal dlffaranclnq on a 
nonunlfora qrld. Tha difference between tha 
achaaaa Ilea In tha aaact detalla of tha 


apatlal dlffaranclnq. Tha flrat achaaa la, 

(v T <r )*|* r i-|’ M t < 7 r* 4 « , )°r‘ ’ 0 <J ” 


with 
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( 22 ) 


Thla la only flrat ordar accurate on non- 
unlfora qrlda. Tha aacond achaaa la aacond 
ordar accurate. 


( r I J ( r ) J 


’t T 


A f 


Tf 
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( 23 ) 


■oth aiaaplaa uaa a eoapu ta t 1 ona 1 doaaln 
with 200 nodaa, and non-unlfora atratchlnq 
auch that r varlaa aa mown In flqura 1. Tha 
Initial condition la a wavapackat located at 



tha etntir of the dooiin, vlth a vavalanqth 
of ipproniMttiv <2 nod* point*. 
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Piquro 4. wavnpackat path In *-i pnaao piano 


Plqura i. Plot of r(kl for both aaaaipl**. 

Plquroo 2 and 1 ahov tho davalopaant of 
tho aoiutlon for tho two aaaapla*. Qualita- 
tively tho behavior in both caooa la tno 
■aao. Por tho flrot 100 ltorationa tho 
vavapacket travola rlqnt. with ineroaoinq 
eoaputat lona 1 vavonuabar aa tha local aoah 
apacinq aa incroaaaa. Aa tha vavonuabar 
incraaaaa paat • / 2 , tha vavapackot ravaraaa 
diractlon, and in tha noat 100 itaratloni 
travala back to tho cantor of tha doaain and 
tho vavonuabar raacnaa a paak of nearly ». 

Tho vavonuabar than baqlna to dacraato vhlla 
tha vavapackot continual aovlnq lift. Attar 
about 300 ltorationa tha vavonuabar la sack 
to t/2, and tha vavapackot ravaraaa diractlon 
aqaln and travala rlqnt back to tha cantor of 
tha doaain, at vhlch tlao It haa tho aaaa 
-vavonuabar It had orlqlnallv. 

Tha analyala of tha flrat aaaapla baqlna 
by caleulatlnq tha dlap.raion ralatton. 

a - -21 a in ( Q/2 ) »lr coa ( a / 2 I a 1 n ( a ) (24) 

O 

ao tha dlaporalon ralation la, 

2 t*n(a/2) • r al.i(a) (25) 

Por a particular value of a tha dlapar- 
alon ralation ahove that thara ara tvo value* 
of a for each value of r. If r la laaa than 
a critical valua r er 1 1 • 2 tan ( a/2 ) , both root* 
ara real vlth ona In tha Interval (0,t/2| and 
tha other In tha Interval (v/2,t). If r la 
qroator tffan r erlt , both root* ara coaplaa. 

A full dlacuaalon of thoao faaturaa la qlvan 
In Ml. Por tho aaaapla praaantad hare 
r<r crlt In tha alddle *0t of tha doaain. 
Plqura 4 ahova tho roal roota a in thla part 
of tha doaain. 

Tha vavapackot aquation* for tho flrat 


aaaapla era. 


11 . 
dn 

r i - 

r coa (a 1 coa 1 ( a/2 ) 

( 24 ) 

11 . 
dn 

dr 

d) 

aln ( a ) coa 1 (Q/2 1 

( 27 ) 

d t 

— • 
dn 

0. 


( 2*) 


equation* (24) and (271 daacrlb* an anti- 
cloekvlaa notion around tho ■-« curve in 
fiqur* 4. At tha turning point* r a r er | t , 
*>«/2 and tha qroup velocity r, la taro. 


contlnuoo aovlnq round tho eurva. 


Plqura 5 ahova a eoaparlaon of tha peti- 
tion and anarqy of tha vavapackot obtained 
froa tha nuaarleal aiparlaant and froa lntaq- 
ratlnq tha vavapackot aquation*. The aqraa- 
aant 1* aurprl ilnqly good conaldarlnq that 
the theory la aoyaptotlc, not aaact. In fact 
tha aayaptotic theory la not atrlctly valid 
at the turninq point*, but follovlnq the 
procedure uaad by Llqnthlll (3) to analyra 
cauatlea (turninq polnta In analytic equa- 
tion*) It la poaaibl* to conatruct a local 
analyala In tho neighborhood of tha turninq 
point, vhlch raaolva* thla difficulty. 

In tha aacond oaaapla, 

a • -21 aln(a/2) ♦ lr coa(Q/2)aln(*) 
o 

♦ co* ( 2 / 2 ) ( coa ( a ) - 1 ) * 0(L-1) (25) 

o 3 

To flrat order In L thl* la tha aaaa aa 
In the flrat aaaapla, ao tha dlaporalon 
aquation and tha aquation* f or cnangai In a 
and a ara all tha aaaa. Tha anarqy aquation 
hovavar la different. 

7 ^ • 2 - 7^1 ' -coa ( # )) coa* (a/2 ( C (30 1 

d n d 3 

Plqura 4 ahova the eoaparlaon between 
the nuaarleal aiparlaant and tha vavapackot 
theory. Tha theory accurately predict* tha 
larqa chanqa in the anarqy of tha vavapackot 
durlnq the oaclllation. 

It light b* argued that to calculate tha 
atablllty of tha finite difference achaaa on* 
ahould raally analyte tha alqanaodaa of the 
finite difference aquation*. In (I) It vaa 
ahovn that for a particular elaaa of prooiaaa 
tho avaraqa decay rat* calculated uainq tho 
vavapackot analyala la aayap t o t 1 ca 1 ly equal 
to the decay rat* of tho alqanaodaa. Uainq a 
alailar analyala It ean bo provad that the 
aaaa la true for thla problaa. and ao alnc* 
tha vavapaciata have taro avaraqa decay rata, 
the alqanaodaa alao have taro decay rat*. To 
daeonetrat* thla thla nuaerlcally, and 
aaaaina tha tranaltlon froa a vavapacket fora 
Into an alqmaoda roproaanta tlon aaaapla 1 









( 3)901 


9 « r> 


NUMERICAL experiment 
wrtvEPACKET Tme:ht 
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NUMERICAL experiment 


• >- • r r r ’ r . . r " - * 

tut’ "« , 't I . . 



Plaura 1. Compariaon between tneory and aapannent naure *• Comparison between theory and eaperioent 

for taaaiple l. tor •*•**>■• 


f 


aa continued far 10,000 Itiratlom. Plaura 
ehewa the aalutlon at warloue ataaaa in Ita 
eve lepaan t . >000 ltaratlona rapraaant 2.1 

erloda of tha aaclllatlon. to in flqura Na) 
ha vavepacket la alternately traveling laft 
hah rlaht. Cradually tha waveoacaet uocoaoa 
tretenad with tha laraaat aaplituda at tha 
!rant af tha wawapackat and a ataadllv jrow- 
«a ‘tall* behind. Thia la a aara aatraaa 
laaapla of a phaneaaaen noted and dlaeuaaad 
ip viehnawetakr [71. In tha original vave- 
laekat tha aaplituda •adulation correaponda 


a aaall perturbation in tha frequency, and 
naequently thara la enerqy aaaorlatad with 
lqhtly hlqhar fraquanclaa which trawaia at 
eaaller freup velocity than tha oejerlty af 
a anarqv. After about 10,000 ltaratlona 
a ana r q y la apraad ihrouqnout tha ration In 
ich r»r er i t . At thla ataqa In tha develop* 
nt tha aalutlon la boat canaidarad ta be a 
a of elqennedea of tha ayataa. Since thara 
a prababiv aavaral elqanaedae with frequen- 
aa claaa ta tha fraquancv of the original 
vapacnat, thara la conildarabla ’baatlnq* 



Continuation of X. 


N- 10000 


„ , . u ill ' * «*IUvA 

- !/ * i I , 'f V V V ; V y IY ;; ' A 



0. c . 50 y 0. '5 1.3 


Flour* 'b. Continuation of isaacl* 1. 
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Abstract 

The initial boundary value problem describing the evolution of unsteady 
linearized perturbations of a steady, unifona subsonic flow is analyzed. The 
eigennodes and eigenf requeni ces of the system are derived and several examples 
are presented to illustrate the effect of different boundary conditions on the 
exponential decay rate of the eigenmodes. The resultant Implications for the 
stability and convergence rates of finite difference computations are 
d. scussed. 
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INTRODUCTION 


In finite difference calculations of steady-state subsonic solutions of 
quasi-one-dimensional and two-dimensional Euler equations using time marching 
methods, it is often observed that when the solution has almost converged to 
steady-state the remaining residual is due to the propagation of low frequency 
waves up and down the domain. These waves are largely unaffected by numerical 
viscosity and are dissipatod hrough the interaction with the inflow and 
outflow boundary conditions. The purpose of this paper is to examine this 
process by analyzing the un~teady linearized perturbations of a one- 
dimensional, steady, uni fo m, subsonic flow. For this linear problem with 
constant coefficients is possible to derive the exact eigenmodes and 
eigenf requencies of the ir'.tial boundary value problem. This is the classical 
technique used to analyze physical and acoustical vibrations in a finite 
domain [5] and more recently used in numerical analysis to examine the 
P-stability of finite difference approximations to scalar equations [2,3’. 
The exponential decay rate of the physical eigenmodes is computed for several 
different sets of boundary conditions commonly used in finite difference 
calculations and the implications for the stability and convergence rates cf 
these calculations are discussed. 

The wellposedness of both the initial boundary value problem (i.b.v.p.) 
and the steady-state boundary value problem (b.v.p.) is discussed briefly. 
The definitive analysis of the i.b.v.p. for multi-dimensional hyperbolic 
systems is given by Kreiss in [4], Olieer and Sundstrom [7], use an energy 
method to establish sufficient conditions for the wellposedness of the Euler 
i.b.v.p. Finally, the wellposedness of the steady-state solution to the 
nonlinear quasi-one-dimensional and two-dimensional Euler equations will be 
discussed in a forthcoming paper by Wornom and Hafez [9]. 
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2. ANALYSIS 

The equation for the unsteady linearized perturbation of a steady, uniform 
one-dlmenslonal flow Is, 



where o, u, p are the perturbation density, velocity and pressure and 
c, u, p are the steady, uniform values. 

The analysis Is greatly simplified by defining che following non- 
di-ensional variables 


p/o 

(2) 

u/c 

(3) 

~ 0 

p/p c~ 

(4) 

X/L 

(5) 

Tc/L, 

(6) 


where 
d ora in 
at x ■ 
The 


- - l h 

c = [ 7 p/p] L is the speed of sound. L 
considered, so in the non-dimensional 
0 and the outflow is at x * 1. 
resultant non-dimensional equation is 


is the physical length of the 
domain the subsonic inflow is 




- 0 , 


(7) 


P 


t 


x 


where 


3 


A - 



1 

M 

1 



I 


and M is the Mach number of the unperturbed flow. 
Equation (7) has wave-like solutions 


provided 


/°\ 

I u I ■ exp i(kx - cot)]l T 
(kA - oj! )U - 0, 


( 8 ) 


(9) 


( 10 ) 


so oj/k is an eigenvalue of A and U is the corresponding eigenvector. 
The three eigenvalues of A and their corresponding eigenvectors are 


X 


1 


M 


* 2 - M + 1 


X 


3 


- M - 1 



(lla,b) 


(12a, b) 


(13a, b) 


A general eigenmode of the initial boundary value problem can be written 
as a sum of the three elgenwaves, 


—iii) t r 


V 


i(y/^)x 


U, 4- a. 


i (■) / ' 2 ) A 


U„ + a. 


i(u/X 3 )x 


(14) 
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The eigenfrequency w and che values of the constants , c^, are 
determined by the three boundary conditions. 

At the Inflow boundary at x ■ 0 there are two boundary conditions which 
when linearized and non-dlnenslonal lzed have the form, 


C 


in 


U 


0 , 


(15) 


where C^ n Is a 2 * 3 matrix. Substitution of (14) into (15) yields the 
equation, 


where 



(16) 


(17) 


A necessary condition for the initial boundary value problem to be 
wellposed Is that the 2*2 matrix 



is nonsingular and so can be Inverted 
the incoming characteristics, as a func 
characteristic. 



to obtain and a^, the values of 

:ion of the value of the outgoing 


1 


5 


Similarly the outflow boundary condition yields one equation of the form 


C U - 0 
out 


(18) 


and substitution of (14) produces 


where 





(19) 


(b 


31 


'32 


b 33* 




out 


K'jAO 

e U. 


i(w/X 3 ) 


( 20 ) 


The second necessary condition for the wellposedness of the Initial 
boundary value problem is that b-jj is nonzero so that the value of the 
incoming characteristic can be determined as a function of and 02 the 
values of the outgoing characteristics. 

Equations (16) and (19) can be written jointly as 


B ( ij) 



- 0 . 


( 21 ) 


To obtain a nontrivial eigenmode B( .) must be singular 
<^2 o j ) must be a corresponding null vector, 
frequenices can be calculated from the following determinant 


and the vector 
Thus the eigen 
equation 


det B(.') - 0. 


( 22 ) 



6 


The matrix B can also be used to examine whether the steady-state 
boundary value problem is wellposed. The three requirements for wellposednc3S 
are that a solution exists, Is unique, and small perturbations In the boundary 
data produce small perturbations In the solution. 

The linearized steady-state boundary value problem has a zero solution and 
this solution is unique provided there are no nonzero solutions to 



(23) 


l.e., provided that B(0) Is non-singular. 

A perturbation of the boundary data leads to an equation of the form 


B(0) 




(24) 


which, provided B(0) is nonsingular, can be solved to obtain 

T 

(u^ a ^ a^) which define the characteristic perturbations of the steady- 

state solution. 

Thus the linearized steady-state boundary value problem Is wellposed if, 
and only if, det B(0) Is nonzero, or alternatively the Initial boundary value 
problem does not have a zero eigenf requency . 


3. EXAMPLES 

(a) Entropy, Enthalpy Specified at Inflow, Pressure at Outflow 


The physical boundary conditions are 


7 


X - 0 


p'/p' Y ■ p/P Y 


’ii i u ^ + i£;.iii-2 +I i 
2 p 2 p 


(25a) 

(25b) 


X - L p' - P, 


(25c) 


where p', u', p' are the unsteady physical variables which are a sura of the 
steady-state and unsteady perturbation variables. The corresponding 
linearized non-dimenslonalized equations are 


x - 0 




(26a) 


x - 1 ( 0 


0 



(26b) 


At x ■ 0 substitution of the eigenvector definitions (lib), (12b), and (13b) 
into the elgennode definition (14) yields 


P 

u 

P 


-lut 

- e 



1 

1 

1 



(27) 


Substitution of this 


equation Into (26a) produces the characteristic inflow 


boundary condition 


8 



,o 

u 

p 


-i-t 

■ e 



1 x \ / a i exp ^ lu)/x i^\ 

1 ” l ) [ '2 ex P^ lj/X 2 ) ) ’ 

1 1/ \n 3 exp(lx/\ 3 )/ 


and substitution into (26b) produces the characteristic outflow bou 
condition 


/' 1 

(0 0 1 ) I 0 1 

\0 1 


1 \ exp(iu/\ l )\ 
-1 I | a 2 exp( Lu/Xj) I 
1' \^ 3 exp( itu/Xj)/ 


- (0 exp(i-/\ 2 ) 


exp( 1 j/\ 3 ) ) 



- 0 . 


Together equations (28) and (30) define the matrix B 


r 

B ('■*)) - [ -1 (>-l)(l+M) 

\ 0 exp( i j/Xj ) 


(Y-l)d-M) 
exp( ij /\ 3 ) 


(28) 


(29) 


•diry 


(30) 


( 31 ) 
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The elgenf requencles are given by 


det B ■ (y— 1)[ (1 ~M) exp(l-/X 2 ) - (1+M) exptiwA^)] - 0, 


““> exp 



1+M 
1 - M ’ 


■“> u> • 


1 - M 
— 1 


-1 log 


(^) 


+ 2n^ 


(32) 


(33) 


(34) 


where n is an Integer. 

Thus there Is an infinite set of discrete elgenf requencles . 
to define a decay rate ? 


def _ 




For this example 


3 

n 





1 1 + M> 

1 - M * 


It is useful 
(35) 


(36) 


The amplitude of the elgennode grows, or decays, as exp(-3t), so the 
requirement for all eigenmodes to decay Is > 0 for every n. In this 

example the requirement Is satisfied and so any Initial disturbance at 
t “ 0 will decay exponentially. 


(b) Mass Flux, F.nthalpv Specified at Inflow, Pressure at Outflow 
The physical boundary conditions are 


X - 0 


,)'u' “DU 



Y-l 


~T~ 


-2 

u 



X - L p' - p. 


(37a) 

(37b) 

(37 C ) 


/ 
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Omitting the algebraic details the resultant matrix B is 


( M 1 + M 

-1 (y-DCI+M) 

0 expd'^Aj) 


(y-1)(1-m) 
exp(ljj/X j) 


( 38 ) 


The elgenf requencies are 


u 

n 



log( 


(1 + M)fl 4 M(y - l) 
(1 - M)[l - M(y - 1) 


|| + 2(n ♦ V 2 


(39) 


The decay rates are 


3 

n 



log 


/(I M) f t + M(y - l)]\ 
M)-[l - M(y " 1)]/ * 


(40) 


(c) Density, Pressure Specified at Inflow, Pressure at Outflow 


The physical boundary conditions are 


X - 0 


X - L 




p - p 


P ■ P 


(41a) 

(41b) 

(41c) 


The matrix B is 


B - 


1 1 

1 1 

expfl./X^) exp(lu)/X^) 


(42) 
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The elgenfrequenlces are 

u> ■ (1 - M^)nn , 
n 


( 43 ) 


and Che decay races are zero. 

Since one of Che elgenfrequenlces is zero Che steady-state boundary value 
problem Is ill-posed, as discussed earlier. 

(d) Denslcy, Veloclcy Specified ac Inflow, Pressure ac Oucflow 
The physical boundary condidons are 

(44a) 
(44b) 
(44c) 


The oacrix B is 

/' 1 \ 

B - I 0 1 -1 j . 

\0 expd'^/Xj) expd-'/X^)/ 

The elgenfrequenlces are 

u n • (1 - M“) (n + V 2 )i 


(45) 


(46) 



and che decay races are -,ero. 

In chls example che steady-scace boundary value problen Is wellposed buc 
because of Che zero dec^y r?tes unsteady oscillations will continue 
Indefinitely without exponential growth or decay. 
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(e) Non-reflecting Boundary Conditions 

The full nonlinear non-reflecting boundary conditions specify entropy 
and the appropriate Rlemann Invariant at the Inflow, and the other Rlemann 
Invariant at the outflow [4] 


The matrix B is 


X = 0 


X - L 


p'/p' Y ■ p/p Y 


. , 2 , _ 2 - 
u + — r c * u + — r c 
7-1 Y-l 


2 , - 2 - 

U - r C * U r C 

Y-l Y-l 


(47a) 

(47b) 

(47 C ) 


-1 

1 

Y-i 


0 

2 


0 

0 


— expdj/Xj) 0 -2 exp^u/A^) 


(48) 


Det B * 0 leads to o = + * which reflects the fact that with these 
boundary conditions the unsteady perturbations become zero after the finite 
tine it ta’Nes for all three characteristic waves to cross the domain once. 


CONCLUSIONS 

The calculation of the exponential decay rates of physical eigenmodes has 
implications for the stability and convergence rates of time-marching finite 
difference computations. If the analytic problem has exponentially increasing 
eigenmodes then for sufficiently fine grid resolution a tine-accurate 
numerical solution will exhibit corresponding exponentially increasing 
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eigenmodes. In a forthcoming paper, Trefethen [8] will prove that for a 
linear constant coefficient system such as this the three conditions: 

(i) Exponentially decaying physical eigenmodes, 

(ii) Dissipative interior numerical scheme, 

(iii) GKS-stable numerical boundary conditions, 
are sufficient to ensure the P-stabllity of a time-marching method for a 
sufficiently fine grid. P-stability was defined by Beam, Warming and Yee [2] 
and corresponds to GKS-stability with the additional requirement that none of 
the numerical eigenmodes increases exponentially. The precise definition of 
the theorem and its proof are given in [8], but In essence the argument is 
that condition (i) ensures that low frequency physical waves decay, while 

conditions (ii) and (iii) ensure the decay of high frequency waves, both 
physi cal and non-physical. 

The exponential decay rates for the physical eigenmodes also provide a 

useful lower limit on the spectral radius of the finite difference time- 

marching procedure. If a physical eigenmode decays as exp(-ot) with 

0 > 0, then for a sufficiently fine grid the corresponding numerical 

eigenmode decays approximately as exp(— 7 r. At) where n is the Iteration 

number and At is the time-step. As the grid is refined with At/Ax held 

constant, At * 0 and so the spectral radius is no less than 
2 

1 - cAt + 0(At ). If o*0, as in example (d), the physical eigenmodes are 
neutrally stable and so the numerical convergence rate towards steady-state is 
due solely to numerical dissipation. If this dissipation is of n L order then 
the corresponding spectral radius is 1 - 0(At n+ ^). Non-reflecting boundary 
conditions as in example (e) clearly give a much faster rate of convergence, 
but in two or three dimensions perfectly non-reflecting boundary conditions do 
not exist and in general the best that can be achieved is that there is zero 
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reflection for locally plane waves propagating in a particular chosen 
direction [1]. 

It is not clear to what extent the conclusions for this model problem, 
with linearized perturbations and constant coefficients, are valid for more 
general flows such as transonic quasi-one-dimensional and two-dimensional 
flows. Nonlinear mechanisms at sonic lines and shocks are undoubtedly very 
important. However the decay to steady-state of low frequency waves will 
still depend on the physical boundary conditions and so this analysis should 
provide insight into the effect of the boundary conditions. 
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Abstract 


The unsteady quasi-one-dimensional Euler equations are solved 
using a conservative box method which is second order accurate and 
requires no non-physical boundary conditions. No artificial viscosity 
is used and so the shock cells and sonic cells require special treatment 
which is related to the behavior of the analytic characteristics. 

Results are given for a converging-diverging channel with a moving shock 
due to the periodic oscillation of the inlet boundary conditions. 
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Notation 

Variables 

A Cross-sectional area of duct 
c Speed of sound 
E Enthalpy 
F Flux vector 
p Pressure 
P Pressure vector 
u Velocity 

U Conservative state vector 
p Density 

Y Ratio of specific heats 
Subscripts 

j Spatial index for discrete variables 
0 Stagnation quantity 
e Exit quantity 

Superscript 

n Iteration number for discrete variables 


& 


C 




1 . Introduction 


The box method was first proposed by Keller [1] for solving 

parabolic equations in which the second order p.d.e. is first rewritten 

as a coupled system of first order p.d.e. 's. It is now widely used for 

solving the boundary layer equations. It also has several attractive 

features for solving hyperbolic systems. It is second order accurate, 

requires no non-physical boundary conditions (such as extrapolation of 

characteristics as required by 3-point differencing schemes ) and for 

the model problem u +cu - 0 it gives the exact answer if cAt/Ax-1. 

t x 

The box method was first used to solve the 1-D Euler equations by 
S. Wornom [2]. Since Wornom was interested in steady-state solutions he 
used a Backward Euler version of the box method and assumed constant 
stagnation enthalpy. For large At this converges rapidly to the steady 
state solution. In supersonic regions artificial compressibilty was 
introduced (as proposed by Eberle [4]) to achieve shock capturing and 
prevent expansion shocks near the sonic point. As a consequence one 
non-physical boundary condition was required at the supersonic outlet 
and this produced a boundary-layer type behavior at the outlet. Wornom 
has also used the time-accurate box method to solve the 1-D unsteady 
Euler equations [3]. In this paper he chose not to use artificial 
compressibilty and consequently did not require any non-physical bound- 
ary conditions but did require special treatment of the sonic cell and 
was unable to handle the shock cell. 

In section 2 of this report the conservative equations are 
derived for subsonic cells (cells for which both the inflow and the 
outflow are subsonic), and supersonic cells (cells for which both the 
inflow and the outflow are supersonic). These equations are identical 
to those used by Wornom [3]. Section 3 derives the physical boundary 
conditions and their numerical implementation. Section 4 discusses the 
difficulties with sonic cells (cells for which the inflow is subsonic 
but the outflow is supersonic), and shock cells (cells for which the 


inflow ia supersonic but the outflow la aubsonic), and their relation to 
the behavior of the analytic charcteriatica. For the 8onic cell the 
difficultiea are reaolved by imposing an additional characteriatic 
equation. For the shock cell a natural form of ahock fitting ia derived 
with the ahock poaition being an additional variable. Section 5 
preaenta a computational example of a convergent-divergent channel with 
a ahock which oacillatea due to a periodic oacillation in the inlet 
atagnation presaure. Finally aection 6 discusses the results and the 
difficulties and current achievements in extensions to the 
two-dimensional Euler equations. 
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2. Equations for Subionic and Supersonic Calls 

The unsteady quasi-or>e-dimensional Eular aquations for a 
variabla araa nozzla ara, 
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An intagral form of Eq. (1) for any computational call is, 
x A . x 

f 3*1 JT| 

/ AU dx ♦ Af| ♦ / At^ dx - 0 (5) 

x j 


This equation remains valid even whan thare is a shock in tha 

ip 

interval [x ,x ] (provided — is correctly represented by a Dirac 
j j+i ax 

delta function) and is tha basis of tha finite difference equations. 


As illustrated in figure 1 the approximations which are 
introduced are that A is piecewise linear between x. and x a . ana U is 


piecewise constant between x_. ^ and 
eq.(5) becomes, 
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where k. -A.«-a3(A. -A.). Now 4r<J is spproxlmatsd by — (U^’-u") 

3*a3 3 3 + i 3 at 3 rr at 3 3 

and tha flux tsra F ^ (and similarly P ^ ) is approximstad at tima laval 

. .. , , „n+9 „n SF^.n+l ,.n. . 

n+e by a linearized expansion F^ U j “ V • ® -1 corresponds to 

tha Backward Euler version while 6*1/2 corresponds to the time-accurate 

box scheme used by Wornom in [3]. Thus the finite difference equations 
which are used are, 
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In matrix form the equations are, 
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3. Boundary Conditions 


a) Supersonic Inlet and Outlet 

At supersonic inlets there are three boundary conditions which 
together totally specify the flow variables at the inlet. At supersonic 
outlets there are no physical boundary conditions, and no boundary 
conditions are required for the numerical solution. 

b) Subsonic Inlet 

At subsonic inlets two boundary conditons are required and they 
are chosen to be specified stagnation pressure and density. For 
programming simplicity these are implemented through the following 
equivalent two conditions. 



which when linearized becomes, 


u, ‘ c2 )" Ac" - (Y-1)u l 4 (ou)" + ( Y-1 ) A(pE)" - -p" + P 0 (p"/p 0 ) Y (15) 
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which when linearized becomes, 
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c) Subsonic Outlet 

At subsonic outlets one boundary condition is required which is 
chosen to be specified exit pressure. 
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n+1 

Pj " P e ( 18 ) 


which when linearized becomes, 

uJ jj “ ( Y~1 ) u j A(pU)” + ( y— 1 ) A(pE)" " P a " Pj 


( 19 ) 
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4. Equations for Sonic and Shock Cells 


In a case in which the flow is supersonic in the entire domain, 
the inlet flow is specified and the matrix equation (7) can be solved 
to obtain AU 2 and then the other AlT by marching downstream. In a case 
in which the flow is subsonic in the entire domain, there are 3J 


variables, (the AU at 1<j<J), 3(J-1) cell equations and 3 boundary 
conditions, two at the inlet and one at the outlet. Thus the number of 
equations equals the number of variables. The system can be written as 
a block tridiagonal set of equations with 3x3 blocks and solved by 
standard methods. Jf 8-0.5 the equations reduce in both cases to those 
used by Wornom [3]. 


Solutions for the above two flow problems are straightforward. 

The difficulties arise with shocks and sonic points. Consider the case 
in which the flow is supersonic at the inlet and subsonic at the outlet. 
There are still 3J flow variables and 3(J-1) cell equations, but now 
there are 4 boundary conditions so there is one more equation than 
variable. This can be understood by considering the behavior of the 
analytic characteristics. The characteristic velocities are u,u±c so 
there are four characteristics entering the shock cell, three on the 
supersonic inflow side of the cell and one on the subsonic outflow side. 
For a steady-state problem the "characteristic information" from the u-c 
characteristics entering the cell on each side must match, but in an 
unsteady problem the mismatch determines the shock velocity. This 
suggests that some form of shock fitting is neccessary. The procedure 
chosen was to define one additional variable, x g , the shock position. 

In the shock cell it is assumed that, 
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as illustated in figure 2, and bo defining X by 
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x (t ) - x + A (x -x ) 
an j 3+1 3 
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the conservation equation with spatial discretization is, 
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and the fully discretized finite difference equation is, 
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Note the two terns involving AX°, the movement in the shock 
position. The first comes from the contribution to 

3 + 1 

AU dx 

j 

due to the movement of the shock, and the other is because the pressure 
jump acts across the shock area which changes in size when the shock 
moves . 



Now that the number of variables equals the number of equations 
the system is well-posed and can be solved. If the shock moves upstream 
past the supersonic node, U at the supersonic node is replaced by U at 
the subsonic node. If the shock moves downstream past the subsonic node 
the opposite is done. One feature of this procedure is that it cor- 
rectly calculates the velocity of a uniform moving shock in a constant 
area duct. Another is that in steady-state solutions the error in the 
shock position is 0(ax j ) and so the global second order accuracy of the 
box scheme is preserved by this shock treatment. 
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In the case in which the flow is subsonic at the inlet, becomes 
sonic at a throat and is supersonic at the outlet, there are 3J vari- 
ables, 3 ( J- 1 ) cell equations and 2 boundary conditions. Thus there are 
one too few equations. An explanation of this is that at the sonic line 
two characteristics emerge with characteristic velocity u-c, one travel- 
ling upstream and one travelling downstream. In the Navier-Stokes 
equations the value propagated along these two characteristics would be 
determined by viscous forces in the neighborhood of the sonic point, but 
in the inviscid Euler equations an extra equation is required to set the 
characteristic value at the sonic point. By diagonalizing the Euler 
equations it can be shown [6] that the equation satisfied at the sonic 
point by the characteristic variable J_ with characteristic velocity 
u-c-0 is 

H- = i£ - cu — - 0 (24) 

3 1 - at at 

Thus the numerical condition which is imposed is that AJ_*0 
on the subsonic side of the sonic cell. When expressed in terms of 
conservation variables this becomes, 

U±lu>] n Ap n - A(pu) n + (y- 1) A ( pE ) n - 0 (25) 

l 2 Jj J j D 3 

One other problem was found in actual computations near sonic 
points. A negative shock from subsonic to supersonic flow is a valid 
solution of the steady state Euler equations, and in practice small 
negative shocks often occurred. According to an inviscid isentropic 
analysis using Riemann invariants [5] these negative shocks should be 
unstable to small perturbations and should become expansion fans. 
Therefore the problem was solved by identifying negative shocks at the 
sonic point, and when they occurred smoothing the two points on either 
aide of the shock. The shock then turned into an expansion fan and 
quickly the solution at the throat became an almost linear expansion 
from subsonic to supersonic flow. 
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5. Example 


The test example is the unsteady problem of the flow through a 
choked converging-diverging nozzle with subsonic inflow and outflow, 
constant exit pressure and stagnation entropy and oscillating stagnation 
pressure. The nozzle area was defined by, 


A ( x ) 
A* 


♦ ( x-0 . 5 ) 2 


(26) 


with 0<x<1 . 50 node points were used with x^ defined by, 


.(ill , f _ lYi . [fill. I 

j l 49j ^49 J 2 1 1 4 {{ 49 J 2 


which gives slightly greater resolution near the throat. The 
stagnation pressure was defined by, 

p - — (p +p ) + — (p -p ) cos(-^-^) 

0 2 max min 2 max min 200 


( 27 ) 


( 28 ) 


with, 


P /P 
max exit 


1 .4 


P min /p .xlt * ’- 2 


(29,30) 


For reference steady flow is choked for Po / ^P ex i t > ^ *235 , and 

figures 3,4 show the steady-state solutions corresponding to p ^ and 

p with the former being unchoked and the latter choked. Figure 5 
max 

shows the unsteady solution which remains choked at all times. The 
throat remains close to sonic throughout the oscillation and so there is 
little oscillation in Mach number between the inlet and the 3hock. The 
shock position varies greatly and so the outflow Mach number oscillates 
considerably. One interesting feature is the "wiggles" in the solution 
near the sonic line. The spatial wavelength of the induced oscillations 
is proportional to the local characteristic velocity. Near the sonic 
line u-c is small, so the wavelength of oscillations of the correspond- 
ing characteristic is small. 
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This problem has a time-step stability limit of approximately, 
u a At/&x<0.2 where u s is the shock speed. The instability occurs when 
the shock is moving downstream. On the subsonic side of the shock cell 


U is defined to be constant n x q <x<x J , . 

3 3+J 

4-1 


In the worst case x„- x J and 

8 j 


so if Ax g > (x^ + ^ -x^ )/2 tii 64% ^ is constant on a negative length. This 

justifies a stabi imit of u g At/Ax<0.5 . i n practice the block 

inversion in the sufc ic solver becomes nearly singular at a lower At. 
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6. Conclusions 


The most important conclusion is that it is possible to solve 
the Euler equations using a truly inviscid numerical method with no 
artificial viscosity and no spurious numerical boundary conditions. 
However this requires shock tracking and special treatment at sonic 
points to replace viscous mechanisms with appropriate inviscid 
conditions. In 1-0 it is relatively easy to produce an efficient 
time-accurate scheme, but in 2-D the problems are much greater. In [3] 
Wornom develops a numerical scheme for steady-state two-dimensional 
Euler flow and presents an example of a supersonic shock reflection 
problem. Unfortunately the method is severely limited by a requirement 
what each of the characteristic velocities u,u±c,v,v+c must not change 
in sign in the domain. In [7] Drela and Giles also develop methods for 
steady-state two-dimensional Euler flow using a conservative streamtube 
formulation. In supersonic applications the solution can be marched 
downstream, and accurately captures shocks without the introduction of 
artificial viscosity. In subsonic applications a very efficient 
relaxation procedure, similar to potential solvers, is used. In 
transonic applications artificial compressibility (very similar to that 
in [2]) is used to capture shocks. Further work is being done to 
improve the performance of the transonic solver. One possibility is a 
special treatment of sonic and shock cells in a manner analogous to that 
used in this paper. However at present no procedure for doing this has 
been found and, looking further ahead to 3-D calculations, it seems 
probable that a shock-capturing scheme will be much easier. 
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Figure 2. Illustration of piecewise constant definition of U in 
a shock cell with shock movement. 
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Figure 5. Mach number at equal intervals during oscillation of 
inlet stagnation pressure 


